function out=modelnash(r,sigma,theta, A)
	size=length(r);
	solq=zeros(size,1);
	u=@(x) 74.1752*x.^(0.6);
	c=@(x) 8.23249 *x.^(1.51678);
	uprime=@(x) 0.6*74.1752*x.^(0.6-1);
	cprime=@(x) 1.51678*8.23249 *x.^(1.51678-1);
	w=@(x) theta*uprime(x)./(theta*uprime(x)+(1-theta)*cprime(x));
	z=@(x) w(x).*c(x)+(1-w(x)).*u(x);
	for j=1:size
		obj=@(x) -(-r(j)*z(x)+sigma*(u(x)-z(x)));
		solq(j)=fminbnd(obj,0.00001,5);
	end
	balances=z(solq);
	out=balances./(sigma*balances+A);
end